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SOLUTION  TO  THE  COMPRESSIBLE  NAVIER-STOKES  EQUATIONS 
OF  MOTION  BY  CHEBYSHEV  POLYNOMIALS  WITH 
IMPLICIT  TIME  STEPPING 


1.  INTRODUCTION 

Gottlieb,  et.  al.  (reference  1)  applied  pseudospectral  methods  to  the  solution 
of  the  one  dimensional  propagating  shock  wave  problem.  They  utilized  a  low  pass 
spectral  filter  which  they  developed  together  with  a  Shuman  filter,  applied  to  the  flow 
on  either  side  of  the  shock  wave  but  not  across  the  shock  front  itself.  The  shock 
location  was  determined  by  examination  of  the  spectral  coefficients.  Since  then,  the 
present  author  has  developed  new  techniques  for  use  with  pseudospectral  methods 
which  have  greatly  increased  their  utility  in  solving  inviscid  flows  with  single  or  multiple 
discontinuities.  Pseudospectral  methods  have  been  used  by  the  author  to  solve  many 
classes  of  complicated  time  dependent  compressible  flows  using  the  full  Euler  equations 
of  motion  (references  2  through  6).  Results  have  shown  that  flow  discontinuities  such  as 
shock  waves  or  contact  surfaces  are  properly  resolved  as  sharp  discontinuities.  Solutions 
for  transonic  airfoil  flows  at  subcritical  and  supercritical  conditions  (reference  6)  were 
obtained  more  recently  and  proved  that  full  pseudospectral  computational  methods 
could  also  successfully  treat  compressible  inviscid  flows  about  non-planar  geometries. 
With  the.  completion  of  the  airfoil  work,  the  author  has  shown  that  pseodospectral 
computational  methods  are  fully  suitable  for  solving  many  classes  of  inviscid,  time 
dependent,  compressible  flows  using  the  Euler  equations  of  motion. 

The  next  logical  step  is  to  turn  to  the  full  viscous  equations  of  motion  namely,  the 
Navier-Stokes  equations.  Orszag  is  the  most  notable  in  the  field  when  one  considers 
stability  and  transition  analyses  of  incompressible  hydrodynamic  boundary  layer  and 
channel  flows  (reference  7  for  example).  He  was  the  first  to  apply  pseudospectral 
methods,  treating  flow  stability  problems  for  laminar  hydrodynamic  boundary  layers. 
However,  no  work  has  as  yet  been  done  for  compressible  flows.  Even  when  one  considers 
the  solution  to  compressible  external  flows,  pseudospectral  methods  have  not  been  used 
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at  all  due  to  difficulties  in  resolving  discontinuities.  The  author  showed  in  reference  8 
that  the  time  dependent,  compressible,  full  Navier-Stokes  equations  of  motion  could  be 
solved  using  full  pseudospectral  methods.  In  that  work,  a  laminar,  oblique  shock  wave 
boundary  layer  interaction  flow  on  a  flat  plate  in  a  supersonic  free  stream  was  computed. 
The  shock  waves  were  maintained  as  sharp  discontinuities  and  the  separation  zone  was 
obtained.  The  spatial  extent  of  the  separation  zone  as  well  as  the  plateau  pressure 
agreed  well  with  experimental  data.  The  only  drawback  was  the  amount  of  CPU  time 
required  to  obtain  convergence.  With  explicit  time  stepping,  and  a  Courant  number 
of  2.5  based  on  the  minimum  time  step  allowable  in  the  entire  computational  domain, 
_fl-ve  hours  of  CRAY-XMP/12  time  (25,000  time  steps)  were  required. 

To  remedy  this  unacceptably  large  machine  time  requirement,  an  implicit  time 
stepping  procedure  was  incorporated  into  the  code.  This  report  presents  results  ob¬ 
tained  with  the  implicit  time  stepping  version  of  the  code  developed  in  reference  8. 


2.  GOVERNING  EQUATIONS 


The  full  two-dimensional,  time  dependent,  compressible  Navier-Stokes  equations 
of  motion  cast  in  conservation  law  form  are  shown  below. 


dU  dE_ 

dt  +  dx 


dF 

+  ^  =  0’ 


‘  where  U  ,  E  and  F  are  vectors  whose  elements  are: 
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where  <7*  ,c7y  axe  the  normal  stresses, rIy  and  ryi  axe  the  sheax  stresses,  A  and  \x  are 
viscosity  coefficients  (Sutherland’s  relation  is  used  since  the  present  work  deals  with 
laminar  flow)  and  k  is  the  coefficient  of  thermal  conductivity.  The  pressure  is  obtained 
from  the  following 
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x  10 
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e  =  r  ~  ~  +  \p(u2  +u2), 


(7-1)  2' 
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The  physical  flow  variables  axe  non-dimensionalized  in  the  following  manner. 
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The  normal  and  shear  stress  terms  axe  non-dimensionalized  by  the  free  stream 
pressure  head,  (/>i?7j )  .  Subscript  one  denotes  free  stream  properties  upstream  of  the 
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incident  shock  wave.  With  respect  to  non-dimensional  flow  variables,  equation  one 
becomes: 
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dt  +  dx  +  dy 
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with  the  Prandtl  number  Pr  defined  by, 
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(4e) 


This  completes  the  non-dimensionalization  of  the  physical  flow  variables  and  the 
conversion  of  the  Navier  Stokes  equations  to  non-dimensional  physical  flow  variables. 
However,  it  still  remains  to  transform  these  equations  into  a  suitable  computational 
space.  This  is  discussed  below.  Several  coordinate  transformations  are  applied  to 
generate  an  appropriate  distribution  of  points  in  the  flow  field.  Appropriate  here  means 
many  points -in  regions  of  large  gradients  and  simultaneously  few  points  in  regions  of 
small  gradients.  The  final  computational  coordinates  are  obtained  using  a  sequence  of 
four  coordinate  transformations.  Namely, 
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The  terms  C\,  C3,  Ax,Ay  and  ,  a  are  transformation  clustering  constants  which 
affect  the  distribution  of  points  in  the  computational  and  physical  domains.  The  final 
form  of  the  Navier  Stokes  equations  becomes, 


Ut  +  Ef  +  Fz  +  H  =  0 


(8) 


where 


U  =  U 


E  =  Elx 


E  EQ^x  +  E^r/y 


H  =  -E&)f  -  [£(<£,)'  +  F&fjy);} 


All  spatial  derivatives  appearing  in  equation  8  axe  calculated  by  pseudospectral 
means.  In  the  present  work  this  involves  the  use  of  chebyshev  polynomials.  The  time 
derivative  Ut  appearing  in  equation  8  is  evaluated  using  finite  differences.  Specifically, 
the  Adams  Bashforth  algorithm  is  used.  The  resultant  difference  form  of  equation  8  is 
given  by, 


*** = * + l6t[di]t  ~  " + ~ 


<9C 


(9) 


The  term  D{j  is  an  artificial  viscosity.  In  the  present  work  fourth  order  artifi¬ 
cial  viscosity  is  utilized  in  x  and  second  order  in  y.  They  are  computed  using  finite 
differences.  The  finite  difference  representations  are  shown  below. 


Di,j  — *  Hx  +  Hy 
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where 


—  —Dx[Ui,j+  2  +  tTi,j- 2  —  4(C^,j+i  +  Uij- 1)  +  6£7jj] 


Py  ~  Dy[Ui+l'j  ~Ui,j  +  U I  — l,j] 


The  terms  Dx  and  Dy  are  smoothing  constants. 


3.  PSEUDOSPECTRAL  METHODS 

Pseudospectral  solution  techniques  involve  the  use  of  series  of  functions  to  represent 
the  global  properties  of  a  flow  field  and  its  spatial  derivatives.  In  the  present  work 
Chebyshev  polynomials  are  used.  They  are  represented  by  Tn(x)  where 


Tn(x)  =  cos  [narccos(x)] 


(11) 


or 


Tn{8)  =  cos  [nd\ 


(12  a) 


where 


8  =  arccos  (x) 


(126) 


A  function  of  a  single  spatial  variable  and  time  such  as  F(x,t)  may  be  represented 
as 


N 

F(x,t)  =  '£An(t)Tn(x) 

n=0 


(13) 


The  time  dependence  is  represented  entirely  in  the  spectral  coefficients  An(t)  while, 
the  spatial  dependence  is  represented  in  the  Chebyshev  polynomials  Tn(x).  The  Cheby¬ 
shev  polynomials  are  evaluated  at  discrete  points  Xj  where 
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(14) 


r  KJ  i 

xi  =  cos  [—  ] 


where  Nx  is  the  total  number  of  modes  used  to  represent  the  spatial  variation  of 
the  function  F(x,t).  The  spatial  derivative  of  the  function  F(x,t)  is  represented  as 


n=0 


where 


Nx 


An(1)(t)  =  pAp(t ) 


p=n+l 


(15a) 


(154) 


and 


p  +  n  =  odd 
Co  =2 
Cn>  o  =  1 

The  An's  are  determined  from  equation  thirteen.  Inverse  FFT’s  are  used  to  obtain 
the  An's  from  the  known  functional  values  F(x,t)  at  the  known  collocation  points  x}. 
The  spectral  coefficients  of  the  spatial  derivative  ,  An{1)  ,  are  determined  from  the 
recurrence  relation  equation  15b.  Direct  FFT’s  are  used  to  evaluate  the  sum  in  equation 
thirteen  to  obtain  the  functional  values  at  t-l-Jt.  The  low  pass  spectral  filter  developed 
by  Gottheb  at  ICASE  is  used  to  damp  spectral  oscillations.  It  is  shown  below. 
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where  K  is  the  spectral  wavenumber  and  Kmax  is  the  maximum  wavenumber 
corresponding  to  the  total  number  of  collocation  points. 


4.  IMPLICIT  TIME  STEPPING  ALGORITHM 

Many  algorithms  exist  to  implicitly  solve  the  Navier-Stokes  equations  of  motion. 
Most  however,  are  not  suitable  for  incorporation  into  a  pseudospectral  code.  The 
reason  is  that  most  rely  on  replacing  various  finite  difference  terms  on  the  right  hand 
side  (ie  the  spatial  discretization  side)  with  their  equivalent  evaluated  at  the  (n+l)’th 
time  step.  Then,  all  the  terms  in  the  finite  difference  representation  of  the  governing 
differential  equation  which  are  evaluated  at  the  (n+l)’th  time  step  are  brought  to 
the  left  hand  side  of  the  equation.  The  resulting  form  of  the  difference  equations  has 
all  (n+l)’th  terms  on  the  left  hand  side  and  all  lower  terms  (ie  n’th,  n-l’st  etc.)  on 
the  right.  The  solution  is  advanced  in  time  by  applying  a  matrix  inversion  to  the 
left  hand  side.  This  procedure  is  not  applicable  herein  because  the  spectral  solver  is 
global  in  nature.  Spatial  derivatives  are  evaluated  row-by-row  or  column- by- column 
(x,y  directions  respectively)  in  one  operation. 

What  is  required  therefore  is  an  implicit  procedure  that  involves  the  time  derivative 
terms  only  since,  they  are  evaluated  using  finite  differences.  The  implicit  procedure 
developed  by  MacCormack  (reference  9)  does  this.  It  is  utilized  herein.  A  brief  outline 
of  MacCormack’s  procedure  will  be  given  below.  The  reader  is  directed  to  reference  9 
for  a  more  detailed  description. 

MacCormack  s  algorithm  is  a  two  step  predictor-corrector  scheme.  Each  step  (pre¬ 
dictor  or  corrector)  involves  the  evaluation  of  a  first  order  accurate  in  time  solution  flux 
or  change,  followed  by  an  implicit  step  which  in  principle  is  stable  for  any  time  step  size 
(  needless  to  say  this  in  practice  means  an  upper  bound  for  the  CFL  condition  of  about 
20  for  laminar  flows).  This  first  flux  calculation  will  be  referred  to  herein  as  the  first 
level  flux  calculation.  In  the  present  work  it  is  calculated  using  the  Adams-Bashforth 
second  order  accurate  finite  difference  time  stepping  algorithm  together  with  the  spec¬ 
tral  evaluation  of  all  spatial  derivatives.  Only  the  predictor  step  is  used  herein  since 
the  time  stepping  algorithm  is  second  order  accurate  and  the  spatial  spectral  algorithm 
is  of  order  accurate  equal  to  the  number  of  modes  used  to  represent  the  flow.  Previous 
work  by  the  author  using  spectral  predictors  and  correctors  has  shown  that  there  is 
essentially  nothing  gamed  by  including  the  corrector  step.  The  CPU  time  spent  on  the 
corrector  i's  wasted  time. 

MacCormack’s  algorithm  is  shown  below. 


(/  -  AiA+^(J  -  A<A+^)«#‘ = % 


(17) 
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Where  the  dot  operator  acts  on  all  terms  to  its  right.  A  and  B  are  matrices  which 
are  themselves  products  of  three  other  matrices.  The  symbol  A denotes  the  forward 
index  finite  difference  operator.  Equation  17  is  solved  in  the  following  manner.  The 
right  hand  side  is  obtained  from  the  Adams-Bashforth  spectral  scheme.  Defining  the 
term  F3  as, 


=  (I  -  At  A+ 


(18) 


where 


f3.  .  =  su^1) 

lyj 

Equation  17  may  be  written  as 


(19) 


\A\ 


(/  "  A<A+AJ^.  =  Ft. 


(20) 


The  procedure  for  solving  equation  20  is  to  use  the  results  of  the  spectral  step,  fj*. 
,and  obtain  F£  .  by  inverting  the  coefficient  matrix  on  the  left  hand  side  of  equation 
20.  Once  the  values  of  F™.  .  are  obtained,  then  a  matrix  inversion  is  applied  to  equation 
18  to  obtain  the  values  of  F3iJ .  Then, the  solution  at  t+dt  is  obtained  as 


(21) 


which  completes  the  procedure  for  one  time  step. 

The  algorithm  includes  the  utilization  of  a  y-direction  artificial  viscosity  term 

whose  magnitude  is  proportional  to  the  solution  fluxes.  It  is  shown  below  in  equation 

22. 


Dy 


(22) 
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This  term  is  required  to  maintain  stability  of  the  implicit  scheme  in  the  early 
stages  of  the  calculation  where  the  fluxes  are  impulsively  large.  As  the  solution  reaches 
convergence,  the  numerator  of  equation  22  reaches  zero  (or  some  small  non-zero  residual 
value)  and  the  dissipation  term  effectively  drops  out. 

The  implicit  part  of  the  overall  procedure  is  the  calculation  of  the  F2  and  F3  vec¬ 
tors.  To  evaluate  these  quantities  at  all  grid  points  would  yield  a  fully  implicit  procedure 
that  would  not  be  the  most  efficient.  For  any  given  value  of  explicit  Courant  number, 
there  result  values  of  explicit  time  step  size  at  each  of  the  grid  points  in  the  computa¬ 
tional  domain.  This  is  due  to  the  variation  of  eigenvalues,  grid  spacing  and  kinematic 
viscosity  coefficient  in  the  computational  field.  These  parameters  change  from  grid 
point  to  grid  point  as  well  as  from  iteration  to  iteration.  The  explicit  procedure  alone 
would  of  course  use  the  global  minimum  value  for  the  entire  field.  The  implicit  proce¬ 
dure  is  implemented  by  specifying  either  an  implicit  courant  number  (greater  than  1.0) 
or  an  implicit  time  step  size  dti  whose  value  is  at  least  larger  than  the  explicit  time  step 
size.  The  explicit  time  step  size  varies  over  the  entire  computational  domain  from  some 
global  minimum  to  some  global  maximum  say,  dti  to  dt2.  The  degree  of  implicitness 
of  the  procedure  is  related  to  the  magnitude  of  dti  relative  to  dt\  and  dt2.  If  dti  is  less 
than  dt\  then  the  implicit  procedure  is  bypassed  and  a  fully  explicit  scheme  results.  If 
dt{  is  greater  than  dt2  then  the  implicit  procedure  is  applied  at  all  points  and  a  fully 
implicit  scheme  results.  In  practice,  dti  is  chosen  to  lie  between  dti  and  dt2.  At  points 
where  dti  is  less  than  dt2  the  implicit  procedure  is  bypassed  while,  at  points  where  dti 
is  greater  than  dti  the  implicit  procedure  is  applied.  Since  the  smallest  field  values 
of  explicit  time  step  size  result  from  the  concentration  of  grid  points  in  the  boundary 
layer,  dti  is  selected  so  that  these  points  are  treated  implicitly  while  others,  which  are 
outside  of  the  boundary  layer,  are  treated  explicitly. 


5.  BOUNDARY  CONDITIONS  AND  TIME  STEP  SIZE 

Both  subsonic  and  supersonic  free  stream  flows  are  treated  in  the  present  work. 
The  type  of  boundary  condition  employed  depends  upon  the  nature  of  the  flow  namely, 
subsonic  or  supersonic.  For  subsonic  flow,  characteristc  boundary  conditions  are  em¬ 
ployed  at  both  inflow  and  outflow  boundaries.  At  subsonic  inflow  boundaries  the  static 
pressure  and  axial  velocity  are  held  fixed  while  at  subsonic  outflow  boundaries  the 
static  pressure  is  held  fixed.  The  remaining  physical  quantities  are  calculated  from 
these  fixed  quantities  together  with  the  values  of  the  characteristics.  Along  the  upper 
computational  boundary,  flow  variables  are  held  fixed  at  free  stream  values  since  the 

boundary  is  placed  far  enough  away  from  the  body  so  that  no  disturbances  propagate 
there. 

_  For  supersonic  flow  problems,  the  inflow  boundary  conditions  are  to  keep  all  flow 
variables  fixed.  Along  the  upper  boundary  flow  variables  were  held  fixed  at  either 
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free  stream  or  post-shock  values  depending  on  the  problem  being  treated.  At  the 
outflow  boundary,  conditions  at  each  point  were  set  to  those  at  the  next  upstream 
point  (zero’th  order  extrapolation)  throughout  the  calculation.  Along  the  bottom  of 
the  computational  boundary  either  a  plane  of  symmetry  or  wall  surface  was  present. 
Reflective  boundary  conditions  were  used  at  points  ahead  of  and  behind  the  body  being 
treated,  namely  =  11^2 •  On  the  body  surface, U{ti  =  — Also,  along  the  entire 
bottom  boundary  —  vi,2  >  Pi,i  —  Pi, 2  and  e,^  =  e,t2,  with  e  denoting  specific 
internal  energy.  This  last  condition  being  applied  along  the  body  surface  since  all  cases 
treated  are  for  an  adiabatic  wall. 

The  time  step  size  was  determined  from  the  following. 


SU  =  [ 


St,  =  [ 


Cn&C  1 

(23a) 

CjvA(2  . 

(23  b) 

irajn 

(23c) 

Where  and  crn  are  the  eigenvalues  in  the  computational  space  at  each  of  the 
grid  points. 

For  the  results  presented  herein,  equation  23  was  utilized  in  the  following  manner. 
First,  a  value  of  explicit  courant  number  was  selected.  Then,  equation  23  was  used  to 
determine  the  resulting  time  step  sizes  at  each  of  the  grid  points.  These  values  were 
then  stored  in  an  array.  As  previously  mentioned  the  implicit  procedure  was  used  by 
either  specifying  an  imphcit  courant  number  and  calculating  the  resulting  implicit  time 
step  size  as  tiie  product  of  the  imphcit  courant  number  and  the  global  minimum  value 
of  explicit  time  step  size  (adjusted  to  reflect  an  explicit  courant  number  of  one)  or 
of  specifying  the  implicit  time  step  size  directly.  In  the  former  method,  the  implicit 
time  step  size  is  allowed  to  vary  from  time  step  to  time  step,  subject  to  the  constraint 
that  the  imphcit  courant  number  is  held  fixed.  While  the  latter  method,  by  fixing  the 
imphcit  time  step  size,  varies  the  imphcit  courant  number  according  to  the  variation  of 
the  eigenvalues  at  the  grid  points.  Note  that  whether  or  not  the  imphcit  procedure  is 
employed,  the  value  of  time  step  size  actually  used  in  equations  17  through  21  at  each 
grid  point  is  the  imphcit  value.  That  value  is  either  a  constant  which  changes  after 
each  iteration  (when  the  imphcit  courant  number  is  specified)  or  is  a  constant  which 
is  held  fixed  throughout  the  entire  computation  (when  the  imphcit  time  step  size  is 
specified  directly). 
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Regardless  of  which  approach  was  utilized,  the  test  used  in  reference  9  was  used 
here  also.  At  each  grid  point,  the  values  of  the  implicit  time  step  size  and  explicit 
time  step  size  for  that  point  are  compared.  For  points  where  the  explicit  time  step 
size  is  greater  than  the  implicit  time  step  size  the  F2  and  Fz  calculations  are  bypassed. 
However,  for  points  where  the  implicit  time  step  size  is  greater  than  the  explicit  time 
step  size  (remember  that  the  explicit  time  step  size  is  different  at  each  grid  point)  the 
F2  and  Fz  calculations  are  performed. 


6.  RESULTS 


Three  types  of  flow  problems  have  been  treated  in  the  present  work.  The  first  is  a 
normal  shock  wave  laminar  boundary  layer  interaction  flow.  The  second  is  a  laminar, 
oblique  shock  wave  boundary  layer  interaction  flow.  The  third  is  the  laminar  flow  over 
a  biconvex  airfoil  at  several  subsonic  free  stream  mach  numbers.  Results  for  each  flow 
will  be  discussed  below. 

The  first  problem  treated  is  a  normal  shock  wave  propagating  into  a  supersonic 
freestream.  The  shock  mach  number  is  4.0  with  respect  to  the  ground.  The  free 
stream  mach  number  is  2.0  with  respect  to  the  ground.  The  body  over  which  the  shock 
propagates  is  a  flat  plate.  The  shock  wave  is  made  stationary  with  respect  to  the  flat 
plate  by  applying  a  velocity  transformation  to  the  entire  flow.  That  transformation 
consists  of  adding  the  negative  value  of  the  shock  propagation  velocity  to  the  shock 
front  and  the  fluid  flow  in  front  of  and  behind  the  shock  wave.  The  resulting  inflow 
is  supersonic  at  a  mach  number  of  2.0  with  the  outflow  also  being  supersonic.  The 
position  of  the  shock  wave  is  chosen  to  be  at  x=0.0.  The  extent  of  the  computational 
domain  is  -1.0  foot  <  x  <  +1.0  foot  and  0.0  foot  <  y  <  0.30  foot.  The  extent  of 
the  flat  plate  is  x=-0.5  foot  to  x=+1.0  foot.  Solutions  for  grid  resolutions  of  64x16 
and  64x32  (x  and  y  directions  respectively)  were  obtained.  Points  were  clustered  in 
the  neighborhood  of  the  plate  surface  to  resolve  the  large  gradients  present.  For  these 
results,  the  implicit  courant  number  was  specified  at  2.0.  The  explicit  courant  number 
was  0.9.  The  free  stream  unit  Reynolds  number  was  purposely  selected  to  be  small  in 
order  to  minimize  the  number  of  grid  points  required  to  adequately  resolve  the  total 

flow  field.  It  was  300,000  per  foot.  This  problem  was  selected  to  serve  as  a  test  to 
check  out  the  implicit  code. 

Results  for  the  64x16  resolution  run  are  shown  in  figures  1  through  7.  These  results 
were  obtained  m  800  iterations.  The  code  required  0.25  cpu  seconds  per  iteration  or 
0.23  x  10  cpu  seconds  per  iteration  per  grid  point.  The  total  cpu  time  required  for 
the  jrun  was  200  seconds  on  the  NRL  CRAY-XMP/24.  Figures  1  through  4  show  full 
field  contour  plots  of  pressure,  velocity  vector  magnitude, mass  flow  per  unit  area  and 
mach  number.  The  sonic  line  is  shown  by  the  dotted  line.  The  computational  grid  used 
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in  this  run  is  shown  in  figure  5  where  constant  coordinate  grid  lines  are  shown.  The 
shock  wave  is  alligned  along  constant  x  coordinate  lines.  In  the  vicinity  of  the  shock 
wave-surface  intercept,  grid  points  are  essentially  evenly  spaced  with  the  point  spacing 
increasing  as  the  inflow  and  outflow  boundaries  are  reached.  The  shock  wave  is  resolved 
sharply  within  two  to  three  cells.  There  is  also  a  small  separation  bubble  present  for 
this  case,  which  is  readily  apparent  in  figures  two  and  three.  The  resulting  surface 
pressure  distribution  ,non-dimensionalized  with  respect  to  the  free  stream  pressure,  is 
shown  in  figure  6.  The  extent  of  the  separation  zone  is  -.125  foot  to  +.125  foot.  The 
sudden  jump  in  pressure  in  the  vicinity  of  the  leading  edge,  at  x=-.5  foot,  is  due  to 
the  lack  of  resolution  of  the  boundary  layer  at  the  leading  edge  and  to  the  change  in 
boundary  condition  of  the  axial  velocity.  The  axial  velocity  goes  from  a  value  of  about 
1.0  just  off  the  leading  edge  to  essentially  zero  at  the  first  point  on  the  flat  plate  surface. 
The  flow  responds  to  this  condition  by  putting  a  weak  shock  or  compression  fan  at  the 
leading  edge  in  order  to  provide  the  physical  mechanism  to  slow  the  flow  down.  This  is 
clearly  evident  in  the  pressure  contour  plot  where  there  are  several  pressure  contours  at 
the  leading  edge  of  the  plate.  A  plot  of  the  residual  time  history  of  the  density  is  shown 
in  figure  7.  The  log  of  the  maximum  density  residual  at  each  iteration  is  plotted  versus 
the  iteration  count.  The  spikes  which  are  present  at  iteration  counts  of  200,  400  and 
600  are  the  result  of  program  restarts  at  these  iterations.  The  code  was  run  800  steps 
hi  increments  of  200.  The  physical  flow  variables  were  stored  with  the  conservative 
variables  being  recalculated  at  each  restart.  Further,  the  implicit  courant  number  of 
2  0  easily  magnifies  any  initial  inaccuracies  which  are  present  at  restart  due  to  the 
recalculation  of  the  conservative  variables.  The  log  of  the  density  residual  has  reached 
-3.5  at  the  800th  time  step.  The  calculation  was  stopped  here  because  there  was  no 
appreciable  change  in  the  surface  pressure  distribution  or  the  extent  of  the  separation 
zone. 

This  same  flow  was  re-run  with  the  grid  resolution  doubled  in  the  y-direction. 
Grid  resolution  of  64x32  versus  64x16  above  was  used.  Results  are  shown  in  figures  8 
through  14.  This  run  required  0.46  cpu  second  per  iteration  or  0.22  x  10“3  cpu  second 
per  iteration  per  point.  Convergence  was  reached  at  3200  iterations  which  required 
14i0  cpu  seconds  or  a  little  over  24  minutes.  The  implicit  courant  number  was  2.0 
as  in  the  case  above.  The  only  difference  between  these  two  cases  is  the  number  of 
chebyshev  modes  used  in  the  y-direction.  At  convergence  the  value  of  the  time  step 
size  for  the  case  above  was  0.60  x  10-3  while  for  this  case  it  was  0.17  x  10~3.  Taking  the 
ratio  of  these  two  and  multiplying  by  800  yields  2823  iterations  required  to  reach  the 
same  integration  time  as  the  above  case  at  the  larger  time  step  size.  The  actual  number 
of  3200  is  very  close  to  this  number  which  takes  into  account  the  effect  of  the  more 
densely  packed  grid  points  in  this  run.  This  is  clearly  evident  if  one  compares  figures 
5  and  12  in  which  the  constant  coordinated  grid  lines  are  plotted  Figure  13  shows 
the  surface  pressure  distribution  for  this  run.  The  extent  of  the  separation  bubble  is 
almost  the  same  as  the  lower  resolution  run.  The  surface  pressure  values  upstream  and 
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downstream  of  the  separation  bubble  are  closer  to  the  respective  pre-  and  post-shock 
values  because  of  the  increased  point  resolution  in  the  y-direction.  Figure  14  shows  the 
residual  time  history.  The  log  of  the  density  residual  has  essentially  plateaued  at  -4.0. 

Based  on  these  two  calculations  it  is  clearly  better  to  have  run  the  coarse  case  first 
then,  after  interpolation,  to  have  continued  the  solution  with  the  higher  grid  resolution. 
Work  using  this  kind  of  approach  as  well  as  multigrid  like  approaches  is  being  presently 
undertaken.  It  was  purposely  not  done  here  because  the  intent  of  this  work  was  to  show 
the  reults  of  incorporating  the  implicit  scheme  into  the  pseudospectral  code  and  the 
utility  of  pseudospectral  computational  methods  developed  in  references  2  through  6 
and  8  for  full,  time  dependent,  compressible,  Navier-Stokes  calculations. 

The  second  case  run  was  a  laminar  oblique  shock  wave  boundary  layer  interaction. 
This  case  was  selected  because  it  had  previously  been  run  with  the  explicit  code  (ref¬ 
erence  8).  Further,  experimental  data  (reference  10)  is  available  for  comparison.  The 
free  stream  Mach  number  is  2.05.  The  free  stream  unit  Reynolds  number  is  695,000 
per  foot.  Sixty  four  chebyshev  modes  are  used  to  model  the  flow  in  the  x  and  y  direc¬ 
tions.  A  six  degree  wedge  is  the  oblique  shock  wave  generator.  The  shock  wave  enters 
the  computational  domain  at  the  left  hand  side  (  supersonic  inflow  )  boundary.  The 
inviscid  reflected  shock  wave  exits  the  right  hand  (supersonic  outflow)  boundary.  The 
upper  computational  boundary  was  positioned  far  enough  above  the  flat  plate  surface 
to  ensure  this.  At  this  grid  resolution  the  implicit  code  required  0.81  cpu  second  per 
iteration  of  CRAY-XMP/24  time.  This  is  equivalent  to  0.192  x  10"3  cpu  second  per 
iteration  per  grid  point.  As  a  comparison,  the  explicit  code  required  0.74  cpu  second 
per  iteration  at  the  same  grid  resolution.  The  explicit  code  required  a  little  over  five 
cpu  hours  of  machine  time  for  convergence  while  the  implicit  code  required  2.7  hours 
with  the  implicit  courant  number  specified  and  2.47  hours  with  the  implicit  time  step 
size  specified.  The  iteration  count  was  12,000  and  11,000  respectively.  From  the  point 
of  view  of  computer  resources  alone,  the  implicit  procedure  indeed  is  a  success,  reduc¬ 
ing  the  computer  time  of  five  hours  to  two  and  a  half  hours.  The  effect  of  the  implicit 
procedure  on  the  solution  will  be  discussed  below. 

Results  for  the  case  where  the  implicit  courant  number  was  specified  are  shown  in 
figures  15  through  40.  Constant  coordinate  grid  lines  are  shown  in  figures  15  through  17. 
The  full  computational  field,  —.2  <  x  <  +.2  foot,  0.0  <■  y  5:  0.3  foot  is  shown  in  figure 
15.  The  clustering  of  points  at  the  plate  surface  is  clearly  evident.  Figures  16  and  17 
show  the  region  of  the  wall  expanded  in  the  y-direction.  Figure  16  shows  0.0  <  y  <  0.08 
foot  and  figure  17  shows  0.0  <  y  <  0.01  foot.  These  ranges  are  used  in  the  other  figures 
and  are  shown  here  so  one  can  visualize  the  number  of  points  that  go  together  with  each 
figure.  The  grid  points  are  essentially  evenly  spaced  in  the  x-direction.  Contour  plots  of 
the  pressure  field  are  shown  in  figures  18  and  19.  Figure  18  is  the  full  field  and  figure  19 
is  a  blowup  of  the  region  0.0  <  y  <  0.08  foot.  The  incident  shock  is  clearly  evident  and 
is  represented  as  a  sharp  discontinuity  even  though  the  computational  coordinates  are 
not  shock  alligned.  The  inviscidly  reflected  shock  wave  has  split  into  two  weak  shocks 
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or  compression  fans.  The  third  arises  from  the  recompression  as  the  flow  re-attaches 
to  the  flat  plate.  Figures  20  through  22  show  details  of  the  separation  bubble  through 
contours  of  velocity  vector  magnitude,  unit  mass  flow  and  mach  number.  The  sonic 
line  is  shown  plotted  as  the  dotted  line.  The  surface  pressure  distribution  is  shown 
in  figure  23.  In  that  figure  the  square  symbols  denote  the  numerical  solution  and  the 
triangle  symbols  the  experimental  data  of  reference  10. 

The  residual  time  history  of  the  solution  is  shown  in  figure  24.  The  high  frequency 
oscillations  which  are  present  over  the  first  two  thousand  time  steps  are  due  to  the  value 
of  implicit  courant  number  chosen.  The  calculation  was  started  with  an  implicit  courant 
number  of  4.0.  This  was  maintained  over  the  first  thousand  time  steps.  The  oscillations 
represent  noise  which  arises  from  the  chebyshev  polynomials  and  does  not  represent 
an  instability.  It  is  typical  of  pseudospectral  chebyshev  solutions  to  display  this  type 
of  behavior  though  of  course  in  practice  the  courant  numbers  used  are  chosen  small 
enough  so  as  to  effectively  remove  these  oscillations  by  minimizing  their  amplitude.  The 
solution  was  kept  because  the  large  values  of  residuals  means  that  solution  changes  are 
propagating  quickly  into  the  field.  The  courant  number  was1  then  reduced  to  three  and 
finally  to  2.0  at  the  start  of  the  1800’th  iteration.  From  the  1800’th  to  the  final  iteration 
the  implicit  courant  number  was  held  at  2.0.  As  mentioned  before,  the  spikes  in  the 
plot  arise  at  solution  restarts.  The  log  of  the  density  residual  reaches  -3.5  at  about 
the  9000’th  iteration  and  rises  slightly  to  about  -3.2  at  the  final  or  12,000’th  iteration. 
The  reason  for  this  rise  is  that  the  dissipation  constants  were  reduced  towards  the 
end  of  the  calculation  to  try  to  minimize  the  effects  of  the  artificial  viscosity  on  the 
solution.  Field  profiles  of  the  physical  variables  at  x=0.025  foot,  roughly  the  middle 
of  the  separation  bubble,  are  shown  in  figures  25  through  40.  There  are  10  points  in 
the  y-direction  spanning  the  separation  bubble  from  the  surface  to  the  bottom  of  the 
separated  shear  layer.  The  solid  vertical  line  in  each  plot  is  the  free  stream  value  of 
the  physical  flow  variable  which  appears  in  the  plot.  This  was  put  in  to  allow  one  to. 
directly  compare  each  of  the  profiles  in  the  viscous  layer  with  the  respective  free  stream 
values.  It  is  seen  from  figures  25  through  27  that  the  profiles  are  smooth  and  that  no 
spectral  oscillations  are  present.  By  looking  at  figures  25  and  27  the  vertical  extent  of 
the  separation  bubble  is  approximately  0.0034  foot.  Note  that  this  is  only  one  percent 
of  the  full  vertical  extent  of  the  computational  boundary.  Further,  the  x  location  of 
these  profiles  is  at  the  maximum  vertical  extent  of  the  separation  zone  so  that  the 
remainder  of  the  separation  bubble  must  be  resolved  in  a  zone  much  less  than  this  one 
percent  value.  The  velocity  profile  which  is  shown  in  figure  26  is  smooth,  with  flow 
reversed  posing  no  particular  computational  difficulty.  Flow  reversal  is  present  from 
the  plate  surface  to  a  value  of  y  just  over  0.0025  foot.  The  extent  of  the  flow  reversal 
zone  is  also  confirmed  by  looking  at  figure  28  where  the  unit  mass  flow  rate  profile  is 
plotted  Figures  29  through  32  show  the  same  profile  plots  over  a  further  extended 
range,  0.0  <  y  <  0.08  foot.  These  figures  provide  a  better  perspective  of  the  size  of 
the  separation  zone  with  respect  to  the  full  field.  Plots  of  these  profiles  over  the  full 
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vertical  extent  of  the  computational  domain  are  shown  in  figures  33  through  36.  No 
oscillations  are  present  at  the  upper  boundary  region.  The  separation  zone  is  no  longer 
distinguishable  in  these  full  field  plots  because  it  is  just  too  small  given  the  scale  of  the 
plots.  However,  the  shear  layer  is  clearly  visible  as  are  the  compression  fans.  Full  field 
(in  the  flow  direction)  profiles  of  density,  velocity  vector  magnitude,  energy  and  unit 
mass  flow  rate  are  shown  in  figures  37  through  40.  The  vertical  extent  of  the  figures  is 
from  the  wall  to  y=0.01  foot.  These  figures  give  a  better  perspective  of  the  numerical 
solution  over  the  full  plate  surface. 

Figures  41  through  63  show  similar  plots  for  the  other  implicit  run.  In  this  run,  the 
implicit  time  step  size  was  held  fixed  throughout  the  computation.  Based  on  an  explicit 
courant  number  of  0.9,  the  resulting  explicit  time  step  size  is  0.22  x  10“4  at  the  start 
of  the  computation  and  0.33  x  10  4  at  the  end.  The  implicit  time  step  size  was  held 
fixed  at  a  value  of  0.80  x  10-4.  This  corresponds  to  implicit  courant  numbers  of  3.21 
and  2.18  respectively.  The  maximum  value  was  about  5.0  early  on  in  the  calculation. 
The  contour  plots  shown  in  figures  41  through  45  are  essentially  the  same  as  those  of 
figures  18  through  22  which  are  results  obtained  with  the  implicit  code  utilized  with 
constant  implicit  courant  number.  Comparing  the  0.0  <  y  <  0.01  foot  profile  plots  of 
the  two  implicit  solutions  (figures  56  through  59,  figures  25  through  28),  the  following 
is  apparent.  The  vertical  extent  of  the  separation  zone  is  the  same.  Surface  values  of 
density  and  energy  are  also  the  same  as  are  the  profiles  in  the  separation  zone.  Velocity 
and  unit  mass  flow  profiles  are  the  same  through  the  separation  zone.  Differences  do 
arise  in  the  shear  layer  which  lies  just  above  the  separation  bubble.  Profiles  of  the  fixed 
time  step  size  run  are  stretched  out  higher  in  the  y  direction  than  the  fixed  courant 
number  profiles.  This  result  is  also  apparent  by  comparing  figures  52  through  55  with 
those  of  figures  29  through  32.  The  maximum  values  of  energy  and  unit  mass  flow 
on  the  profiles  are  however  the  same  for  both  solutions  in  this  vertical  range  off  the 
plate.  Comparing  the  full  field  profiles,  figures  48  through  51  with  figures  33  through 
36  the  vertical  extent  of  the  entire  viscous  layer  (separation  and  shear  zones)  ranges 
through  a  y  value  of  0.10  foot  for  both  implicit  solutions.  This  is  one  third  of  the  full 
vertical  extent  of  the  computational  domain.  The  peak  value  of  energy  reached  in  the 
viscous  layer  is  slightly  larger  for  the  fixed  time  step  size  run.  The  axial  velocity  profiles 
are  essentially  the  same  as  are  the  unit  mass  flow  profiles.  As  previously  mentioned 
differences  do  exist  in  the  density  profiles  in  the  lower  portion  of  the  shear  layer,  just 
above  the  separation  zone.  The  peak  values  of  density  in  the  viscous  laver  are  1.57  for 
the  fixed  courant  number  run  and  1.55  for  the  fixed  time  step  size  run.* 

Some  graphical  results  of  reference  8  (explicit  code)  are  shown  in  figures  64  through 
69  -or  comparison.  The  magnitude  of  the  artificial  smoothing  required  for  stability  was 
much  less  in  the  implicit  code  than  in  the  explicit  code.  The  magnitude  of  the  fourth 
order  smoothing  constant  for  the  x-direction  being  0.03  and  0.002  for  the  explicit  and 
mp™  respectively.  For  the  y-direction  the  explicit  code  smoothing  constant 
was  0.0013  for  a  fourth  order  scheme.  It  was  necessary  to  use  a  second  order  scheme  for 
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the  implicit  code  in  the  y-direction.  The  smoothing  constant  was  0.0007.  The  overall 
effect  of  this  can  be  seen  by  comparing  figures  18,  41  and  64.  The  three  plots  are 
grouped  together  in  Figure  70  for  ease  of  comparison.  These  are  the  full  field  pressure 
contours  of  the  implicit  solutions  and  the  previously  obtained  explicit  solution.  The 
shocks  in  the  explicit  solution  are  a  little  more  smeared  than  those  of  the  implicit 
solutions.  This  effect  is  also  evident  in  the  separation  bubble  whose  surface  extent 
is  somewhat  larger  in  the  explicit  calculation  than  in  the  implicit  calculations.  The 
resulting  surface  pressure  distributions  are  grouped  together  for  ease  of  comparison  in 
figure  71.  In  figure  71  the  square  symbols  are  the  numerical  solution  results  and  the 
triangle  symbols  are  the  experimental  data  of  reference  10.  The  mach  number  contour 
plots  are  grouped  together  in  figure  72. 

The  order  of  the  three  surface  pressure  plots  in  figure  71  is  the  fixed  courant  number 
explicit  code  run  (top  figure),  fixed  implicit  courant  number  implicit  code  run  (middle 
figure)  and  the  fixed  time  step  size  implicit  code  run  (  bottom  figure).  When  comparing 
the  plateau  pressures  the  best  results  are  obtained  with  the  explicit  code.  The  explicit 
code  gives  a  value  of  about  1.28  as  compared  to  the  experimental  value  of  1.25.  The 
fixed  courant  number  implicit  run  gives  1.37  while  the  fixed  time  step  size  implicit  run 
gives  1.34.  In  terms  of  percentages  this  becomes  2.4, 9. 6  and  7.2  percent  respectively 
above  the  experimentally  measured  value.  In  terms  of  the  extent  of  the  separation 
bubble  on  the  plate  surface,  all  three  give  essentially  the  same  recompression  point 
locations.  However,  while  the  fixed  courant  number  implicit  run  and  the  explicit  code 
run  give  similar  locations  for  the  separation  point  (  x  =  -.05  foot  )  the  fixed  time  step 
size  implicit  run  yields  a  separation  point  of  x=-  0375  foot.  Based  on  the  experimental 
surface  pressure  distribution,  the  separation  point  lies  somewhere  between  -0.05  and 
-0.0375. 

The  mach  number  contour  results  of  the  three  cases  are  shown  in  figure  72.  The 
three  solutions  are  essentially  the  same.  However,  as  discussed  above,  there  are  some 
differences  in  the  size  of  the  separation  bubble.  The  sonic  lines  of  the  implicit  code  runs 
are  very  similar,  with  differences  arising  in  the  region  behind  the  reattachment  point. 
The  fixed  courant  number  case  line  at  a  position  higher  off  the  plate  surface  than  the 
fixed  time  step  size  implicit  run.  The  fixed  time  step  size  result  more  closely  matches 
the  explicit  code  result. 

The  third  and  final  flowfield  computed  in  the  present  work  is  the  subcritical  and 
supercritical  flow  over  a  circular  arc  airfoil  of  five  percent  half  thickness  ratio.  Free 
stream  mach  numbers  were  0.70  and  0.84  respectively.  These  cases  were  previously  run 
using  the  pseudospectral  euler  code  that  the  present  author  reported  on  in  reference 
6.  For  the  euler  solution  the  flow  over  the  airfoil  was  supercritical  at  the  free  stream 
mach  number  of  0.84. 

The  subcritical  case  was  the  first  airfoil  case  to  be  run  using  a  full  pseudospectral 
Navier-Stokes  code,  the  implicit  code  in  this  case.  The  implicit  time  step  size  was  held 
fixed  at  0.25~3  while  the  explicit  time  step  size  (based  on  a  courant  number  of  0.9)  was 
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0.313  x  10~3.  The  resulting  implicit  courant  number  was  about  0.72,  meaning  that  this 
run  was  actually  performed  explicitly.  This  was  purposely  done  to  check  out  the  code  for 
this  first  non  planar  geometry  case.  MacCormack’s  y-direction  dissipation  scheme  was 
employed  in  this  run.  The  grid  resolution  was  128  modes  in  the  x-direction,  clustered 
about  the  leading  and  trailing  edges,  and  32  modes  in  the  y-direction,  clustered  about 
the  airfoil/plane  of  symmetry  surface.  This  run  required  0.88  cpu  seconds  per  iteration 
on  the  NRL  CRAY-XMP/24  computer.  At  the  present  grid  resolution  this  works  out 
to  be  0.207  x  10-3  cpu  second  per  iteration  per  grid  point.  Convergence  was  achieved 
after  5000  time  steps,  or  a  machine  time  of  73  minutes.  The  free  stream  mach  number 
was  0.70  with  the  free  stream  unit  Reynolds  number  specified  at  200,000  per  foot. 
Again,  this  low  value  was  chosen  specifically  to  allow  the  numerical  resolution  of  the 
field  with  as  few  modes  as  possible  in  the  y-direction.  Figures  73  through  75  show  plots 
of  constant  computational  coordinates.  The  clustering  is  readily  apparent.  Figure  73 
shows  the  full  computational  field  which  ranges  in  from  -5.0  chord  lengths  to  +5.0  chord 
lengths.  The  airfoil  is  centered  at  x=0.0  with  the  leading  and  trailing  edges  located 
at  x=-0.5  and  x=+0.5  chords  respectively.  The  vertical  extent  of  the  computational 
boundary  is  two  chord  lengths.  Figures  74  and  75  are  enlargements  of  the  area  off  the 
airfoil/plane  of  symmetry  and  extend  respectively,  0.5  and  0.2  chord  lengths  off  the 
plane  of  symmetry. 

Figures  76  through  79  show  full  field  contours  of  pressure,  velocity  vector  magni¬ 
tude,  unit  mass  flow  rate  and  mach  number.  Figures  80  through  83  show  enlargements 
of  the  area  in  the  vicinity  of  the  airfoil  surface.  Figure  84  shows  the  pressure  contours 
and  surface  pressure  distribution  of  the  euler  solution  for  comparison.  It  is  symmetric 
*  with  respect  to  the  airfoil  midchord  since  the  flow  is  inviscid  and  the  airfoil  is  symmetric 
about  its  midchord.  The  largest  difference  between  the  inviscid  and  viscous  solution 
contours  is  in  the  midchord  to  trailing  edge  region.  Figure  85  shows  the  plot  of  the 
pressure  along  the  bottom  of  the  computational  boundary  namely,  along  the  plane  of 
symmetry.  In  comparing  the  Navier-Stokes  and  Euler  solutions  for  pressure  distribu¬ 
tion,  several  differences  are  clear.  The  leading  edge  pressure  coefficient  is  much  larger 
for  the  viscous  solution,  1.35  versus  0.55  respectively.  The  upstream  range  of  influence 
of  the  airfoil  leading  edge  is  about  the  same  for  both  solutions.  The  pressure  coefficient 
has  reached  zero  at  about  1.5  chord  lengths  in  front  of  the  leading  edge.  The  minimum 
pressure  coefficient  on  the  airfoil  surface  occurrs  at  midchord  and  is  -0.55.  While,  for 
the  Navier-Stokes  solution  it  occurrs  at  the  trailing  edge  and  is  approximately  -0.45. 
On  the  downstream  side  of  the  airfoil  the  inviscid  solution  decays  to  free  stream  within 
about  a  half  chord  downstream  of  the  trailing  edge.  The  viscous  solution  disturbances 
extend  to  about  three  chord  lengths  downstream.  Figures  86  an  87  show  the  plane  of 
symmetry  mach  number  and  density  distributions.  The  dotted  line  in  figure  86  is  at  a 
value  of  mach  number  equal  to  the  free  stream  value  of  0.70.  At  1.5  chord  lengths  ahead 
of  the  leading  edge  the  flow  has  returned  to  free  stream.  The  mach  number  attains  a 
value  of  0.72  (the  density  floats  and  is  calculated  from  characteristics),  slightly  higher 


19 


than  the  true  free  stream  value  of  0.70.  Along  the  airfoil  surface  the  mach  number 
is  essentially  zero  since  the  numerical  values  of  the  flow  velocities  are  essentially  zero. 
(Typical  values  are  1.0  x  10-4,  and  1.0  x  10-6  for  u  and  v  respectively.)  The  down¬ 
stream  mach  number  reaches  0.65  at  the  outflow  boundary.  That  this  is  a  lower  value 
than  free  stream  is  to  be  expected  since  the  mach  number  is  a  measure  of  the  total 
pressure  in  the  flow  (at  constant  static  pressure).  The  viscous  solution  must  include 
a  loss  in  total  pressure  to  correspond  to  the  presence  of  viscous  drag.  The  surface 
density  distribution  is  qualitatively  similar  to  the  surface  pressure  distribution.  The 
slight  oscillation  at  the  trailing  edge  is  due  to  the  change  in  surface  boundary  condition 
that  occurrs  in  going  from  the  last  point  on  the  airfoil  surface  to  the  first  point  off  the 
airfoil  surface.  Expanded  plots  of  the  above,  in  the  range  -1.0  <  x  <  +1.0  chord  are 
shown  in  figures  88  through  90.  Full  vertical  field  plots  of  velocity,  energy  and  unit 
mass  flow  rate  are  shown  in  figures  91  through  93.  The  solution  is  seen  to  be  smooth 
with  the  effect  of  the  viscous  layer  clearly  evident  at  the  airfoil  surface.  The  residual 
(log  density)  iteration  time  history  for  the  run  is  shown  in  figure  94.  This  run  was 
completed  in  five  steps  of  one  thousand  iterations  each.  The  log  of  the  density  residual 
monotonically  decreases  to  about  -4.3  at  3500  iterations  where  it  remains  flat  through 
the  5000  th  iteration.  Jt  siiould  be  noted  that  the  surface  pressure  distribution  of  the 
airfoil  omy  converged  in  approximately  one  half  of  the  total  iterations,  or  about  2500. 
The  remainder  of  the  work  was  required  to  obtain  convergence  in  the  field. 

The  subcritical  run  just  discussed,  which  was  effectively  run  fully  explicitly  above, 
was  re-run  again  at  fixed  implicit  time  step  size  but  at  an  implicit  courant  number 
greater  than  one.  Once  again  the  explicit  courant  number  was  0.9,  which  yielded  a 
time  step  size  of  0.314  x  10-4  The  implicit  time  step  size  was  maintained  at  0.6  x  10-4 
v  hich  represented  an  implicit  courant  number  of  up  to  2.25  at  the  first  several  hundred 
iterations  and  a  value  of  1.72  for  the  most  part  thereafter.  All  other  conditions  were  just 
as  those  in  the  case  just  discussed.  Full  convergence  was  obtained  in  2500  iterations. 
The  case  took  0.88  cpu  seconds  per  iteration  or  0.207  x  10-3  cpu  seconds  per  point 
per  iteration.  Total  NRL  CRAY-XMP/24  cpu  time  was  37  minutes  as  compared  to 
73  minutes  for  the  implicit  code  used  in  full  explicit  mode.  Figures  95-98  and  99-102 
show  full  field  and  expanded  field  contours  of  pressure, velocity  vector  magnitude,  unit 
mass  flow  and  mach  number.  Figures  103  through  105  show  full  plane  of  symmetry 
(-5.0  <  x  <  +5.0)  and  figures  106  through  108  expanded  plane  of  symmetry  (-1.0  < 
x  <  +1.0)  plots  of  the  pressure,  mach  number  and  density  distributions.  There  are 
clear'y  differences  in  the  field  contours.  There  are  major  differences  in  the  flow  variables 
along  the  airfoil  surface  itself.  The  mach  number  distribution  on  the  airfoil  surface  is 
uniformly  larger  for  the  implicit  run  (except  right  at  the  trailing  edge).  Since  the 
density  distributions  are  similar,  the  majority  of  this  difference  is  due  to  differences 
in  the  velocity  distribution  along  the  surface.  The  implicit  run  yields  larger  surface 
velocities  which  translate  directly  into  lower  surface  pressure  values.  This  can  be  seen 
directly  by  comparing  Figures  88  and  106.  The  surface  pressure  distribution  of  the 
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explicit  code  solution  decreases  uniformly  along  the  airfoil  surface  until  the  trailing 
edge,  thereafter  recovering.  The  implicit  code  solution  has  a  wide  minimum  just  aft  of 
the  midchord  position,  recovering  from  the  three  quarter  chord  point  on  downstream. 
Figures  109  through  112  show  the  full  field  profiles  of  density,  velocity,  energy  and 
unit  mass  flow  rate.  The  log  density  residual  time  history  is  shown  in  Figure  113. 
The  log  of  the  residual  decreases  monotonically  to  a  value  of  -4.0  (slightly  larger  than 
the  -4.3  value  for  the  case  where  the  implicit  code  was  run  in  explicit  mode).  The 
artificial  smoothing  constants  were  kept  the  same  in  both  runs.  The  difference  in  the 
two  solutions  is  apparently  due  to  the  differing  nature  of  the  schemes  themselves  as 
regards  dissipative  characteristics  and  the  like. 

The  final  case  considered  was  a  full  Navier-Stokes  rerun  of  what,  for  the  euler 
equations,  was  a  supercritical  case.  All  grid  parameters,  airfoil  geometry  etc.  were 
kept  the  same.  The  free  stream  mach  number  was  increased  to  0.84.  MacCormack’s 
y-direction  smoothing  was  utilized.  This  case  was  run  using  fixed  implicit  time  step  size 
(0.4  x  10-3)  at  an  explicit  courant  number  of  0.9  (dt  =0.365  x  10-3)  which  results  in 
an  implicit  courant  number  of  0.985.  The  free  stream  unit  Reynolds  number  was  once 
again  200,000  per  foot.  At  this  low  value  no  shock  wave  was  present  in  the  viscous 
solution.  Full  field  contours  of  pressure,  velocity  vector  magnitude,  unit  mass  flow 
rate  and  mach  number  are  shown  in  figures  114-117.  Figures  118  through  121  show 
expanded  plots  of  the  region  near  the  airfoil.  Figures  122  through  124  show  the  plane  of 
symmetry  distributions  of  pressure,  mach  number  and  density.  Figures  125  through  127 
show  expanded  (-1.0  <  x  <  +1.0)  plots  of  these  distributions.  Full  field  profile  plots  of 
density,  axial  velocity,  energy  and  unit  mass  flow  rate  are  shown  in  figures  128  through 
131.  The  profiles  are  smooth  and  the  viscous  layer  is  clearly  shown.  The  density 
residual  time  history  is  shown  in  Figure  132.  The  log  of  the  residual  monotonically 
decreases  to  a  value  of  -4.4.  It  is  seen  that  the  slope  of  the  curve  is  still  negative  at 
the  point  where  the  calculation  was  stopped.  As  in  the  case  discussed  earlier  the  airfoil 
surface  pressure  has  converged  in  about  half,  in  this  case  2000,  the  number  of  iterations 
required  for  field  convergence.  The  machine  time  to  reach  2000  iterations  is  29  minutes. 
The  euler  solution  surface  pressure  distribution  and  field  pressure  contours  are  shown 
in  figures  132  and  133  for  comparison  with  the  present  Navier-Stokes  solution. 


7  CONCLUSIONS 


Three  classes  of  laminar  viscous  flows  have  been  solved  in  the  present  work. 
Namely,  a  normal  shock  wave  propagating  in  a  supersonic  free  stream,  an  oblique 
shock-laminar  boundary  layer  interaction  and  laminar  airfoil  flow.  This  work  repre¬ 
sents  the  first  time  full  pseudospectral  numerical  methods  have  been  utilized  to  solve 
the  full,  time  dependent,  compressible  Navier-Stokes  equations  of  motion  for  airfoil 
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flows.  The  purpose  of  the  present  work  was  to  develop  a  method  to  accelerate  the 
convergence  of  the  author’s  explicit  (full  pseudospectral)  Navier-Stokes  code,  reference 
8.  To  that  end,  MacCormack’s  implicit  scheme  was  incorporated  into  the  code.  The 
implicit  code  was  checked  out  against  a  generic  laminar  flow  problem,  the  normal  shock 
wave  boundary  layer  interaction  flow,  and  worked  well.  Results  for  this  problem  clearly 
showed  the  effect  of  grid  resolution.  At  64x16  grid  resolution  the  solution  converged 
in  800  time  steps  at  200  cpu  seconds  of  CRAY-XMP/24  machine  time.  At  64x32  grid 
resolution,  3200  time  steps  or  24  minutes  of  machine  time  were  required.  The  difference 
being  solely  due  to  the  extra  constraint  on  the  time  step  size  imposed  by  the  addition 
of  points  in  the  vertical  direction.  The  ratio  of  the  time  step  sizes  actually  used  was 
approximately  equal  to  to  4,  the  ratio  of  the  time  steps  required  for  convergence  for  the 
two  solutions.  This  with  the  respective  time  steps  based  on  an  implicit  courant  num¬ 
ber  of  2.0.  The  implication  of  all  of  this  is  quite  clear.  While  for  this  relatively  simple 
problem  the  advantage  of  the  pseudospectral  approach  (requiring  only  few  grid  points) 
is  clear,  the  overall  point  is  that  it  would  be  more  efficient  for  the  more  complicated 
general  problem  to  have  run  the  solution  at  the  lower  number  of  grid  points  first.  Then, 
at  an  appropriate  point  (density  residual  reaching  some  pre-set  trigger  value  for  exam¬ 
ple)  interpolating  the  solution  to  the  larger  grid  and  resuming.  Multigrid  techniques 
do  this  and  future  work  will  be  directed  in  determining  the  utility  of  such  methods  in 
the  context  of  the  pseudospectral  environment. 

The  second  flow  treated  in  the  present  work  was  the  oblique  shock  wave  laminar 
boundary  layer  interaction  flow  treated  in  Reference  8  with  the  explicit  code.  This  case 
was  re-run  here  to  determine  the  benefits  of  using  the  implicit  code.  All  conditions 
were  kept  the  same  in  order  to  accurately  compare  the  two  codes.  The  explicit  code 
required  25,000  iterations  at  a  courant,  number  of  2.5  based  on  the  global  field  minimum 
time  step  size  and  used  5  hours  of  CRAY-XMP/12  machine  time.  The  implicit  code 
required  11,000  iterations  at  a  courant  number  range  of  2.15  to  3.0  (recall  this  case 
was  performed  with  fixed  implicit  time  step  size)  and  used  2.47  hours  of  machine  time 
on  the  NRL  CRAY-XMP/24.  From  the  point  of  view  of  machine  time,  the  utilization 
of  the  implicit  procedure  developed  by  MacCormack  was  a  complete  success,  cutting 
the  machine  time  to  reach  convergence  in  half.  The  shock  waves  and  compression  fans 
axe  more  sharpie  resolved  with  the  implicit  code  due  to  the  smaller  artificial  damping 
required  to  maintain  numerical  stability.  There  were  however  differences  between  the 
exphcit  and  implicit  solutions.  The  major  one  being  that  the  plateau  pressure  was 
more  overpredicted  with  the  implicit  code  than  with  the  explicit  code  which  essentially 
matched  the  experimentally  measured  surface  pressure  distribution.  At  this  point  it 
appears  ,and  this  is  only  a  supposition,  that  the  inherent  numerical  nature  of  the 
implicit  procedure  is  responsible.  If  this  is  true,  then  the  only  way  to  implement  it 
would  be  to  run  the  implicit  code  to  near  convergence  and  then  complete  the  run  with 
the  implicit  portion  of  the  code  turned  off. 
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The  third  and  final  case  treated  herein  was  the  laminar  circular  arc  airfoil.  The 
subcritical  (free  stream  mach  number  of  0.70)  and  supercritical  (free  stream  mach 
number  of  0.84)  cases  solved  for  were  selected  as  generic  cases  to  determine  how  the 
viscous  pseudospectral  code  would  respond  to  a  solution  of  a  flow  over  a  non  planar 
geometry  with  a  resulting  non  zero  axial  pressure  distribution.  While  the  implicit  code, 
run  in  an  explicit  manner  took  73  minutes  of  machine  time,  the  implicit  run  took  only 
37  minutes.  At  the  low  value  of  free  stream  unit  Reynolds  number  selected  (200,000 
per  foot)  the  euler  supercritical  case  was  no  longer  supercritical.  Of  course  at  a  more 
physically  meaningful  value  of  2  or  3  million  this  would  not  necessarily  be  the  case. 

In  summary  then,  the  implicit  procedure  incorporated  into  the  author’s  explicit 
pseudospectral  Navier-Stokes  code  does  indeed  reduce  significantly  the  machine  time 
required  to  obtain  numerical  convergence.  Questions,  however  remain  as  to  its  actual 
usefulness  since  the  implicit  solutions  differ  from  the  explicit  solutions.  As  a  matter 
of  implementation  the  implicit  procedure  could  indeed  be  shut  off  prior  to  final  con¬ 
vergence.  The  remaining  portion  of  the  calculation  completed  explicitly.  However,  a 
possibly  more  promising  avenue  of  research  lies  in  the  incorporation  of  multigrid  tech¬ 
niques  into  the  pseudospectral  code.  This  offers  a  potentially  much  larger  payoff.  Work 
is  currently  proceeding  along  these  lines. 
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Figure  1  -  Full  Field  Pressure  Contours 
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Figure  3  -  Full  Field  Unit  Mass  Flow  Contours 
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Figure  4  -  Full  Field  Mach  Number  Contours 
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Figure  6  -  Surface  Pressure  Distribution 
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Figure  8  -  Full  Field  Pressure  Contours 
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Figure  9  -  Full  Field  Unit  Mass  Flow  Contours 
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Figure  10  -  Full  Field  Velocity  Vector  Contours 
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Figure  11  -  Full  Field  Mach  Number  Contours 
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Figure  12  -  Constant  Coordinate  lines 
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Figure  13  -  Surface  Pressure  Distribution 
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Figure  15  -  Full  Field  Constant  Coordinate  lines 
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Figure  17  -  Constant  Coordinate  lines 


Figure  18  -  Full  Field  Pressure  Contours 
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Figure  19  -  Pressure  Contours 
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Figure  20  -  Velocity  Vector  Contours 
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Figure  21  -  Unit  Mass  Flow  Contours 
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Figure  22  -  Mach  Number  Contours 
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Figure  23  -  Numerical  versus  Experimental  Pressure  Distribution 
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Figure  24  -  Residual  Iteration  Time  History 
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Figure  25  -  Density  Profile 
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Figure  26  -  Velocity  Profile 
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Figure  27  -  Energy  Profile 


Figure  28  -  Unit  Mass  Flow  Profile 
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Figure  29  -  Density  Profile 
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Figure  30  -  Velocity  Profile 
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Figure  31  -  Energy  Profile 
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Figure  32  -  Unit  Mass  Flow  Profile 
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Figure  33  -  Density  Profile 
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Figure  34  -  Velocity  Profile 
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Figure  35  -  Energy  Profile 
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Figure  36  -  Unit  Mass  Flow  Profile 
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Figure  37  -  Full  Field  Density  Profiles 


U/Ul  VS.  Y 

Ml -2. 05,  RE1-695K/FT,  GRIO:  61  X  61,  1TLOCL-1 


o 


Figure  38  -  Full  Field  Velocity  Profiles 
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Figure  39  -  Full  Field  Energy  Profiles 
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Figure  40  -  Full  Field  Unit  Mass  Flow  Profiles 
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Figure  41  -  Full  Field  Pressure  Contours 
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Figure  42  -  Pressure  Contours 
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Figure  43  -  Velocity  Vector  Contours 
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Figure  44  -  Unit  Mass  Flow  Contours 
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Figure  45  -  Mach  Number  Contours 
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Figure  46  -  Numerical  versus  Experimental  Pressure  Distribution 
(Square  Symbols-Numerical  Results;  Triangle  Symbols-Experiment) 
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Figure  47  -  Residual  Iteration  Time  History 
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Figure  48  -  Density  Profile 


48 


U/Ui  VS.  Y 

XSTfl  -  0.025  FOOT 


Figure  49  -  Velocity  Profile 
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Figure  50  -  Energy  Profile 
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Figure  51  -  Unit  Mass  Flow  Profile 
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Figure  52  -  Density  Profile 
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Figure  53  -  Velocity  Profile 
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Figure  54  -  Energy  Profile 
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Figure  55  -  Unit  Mass  Flow  Profile 
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Figure  56  -  Density  Profile 
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Figure  57  -  Velocity  Profile 
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Figure  58  -  Energy  Profile 
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Figure  59  -  Unit  Mass  Flow  Profile 
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Figure  60  -  Full  Field  Density  Profiles 
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Figure  61  -  Full  Field  Velocity  Profiles 
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Figure  62  -  Full  Field  Energy  Profiles 
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Figure  63  -  Full  Field  Unit  Mass  Flow  Profiles 
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Figure  64  -  Full  Field  Pressure  Contours 


PRESSURE  CONTOURS 

EXPLICIT  CCOCs  M|-2. 05,  RE1-695K/FT,  GRID:  64  X  64 


CD 

O 


Figure  65  -  Pressure  Contours 
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Figure  66  -  Velocity  Vector  Contours 


Figure  67  -  Unit  Mass  Flow  Contours 
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Figure  69  Numerical  versus  Experimental  Pressure  Distribution 
(Square  Symbols-Numerical  Results;  Triangle  Symbols-Experiment) 
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Figure  70  -  Pressure  Contours:  Explicit  versus  Implicit  Results 
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Figure  71  -  Pressure  Distribution:  Explicit  vs.  Implicit  Results 
(Square  Symbols-Numerical  Results;  Triangle  Symbols-Experiment) 
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Figure  72  -  Mach  Number  Contours:  Explicit  vs.  Implicit  Results 
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Figure  73  -  Full  Field  Constant  Coordinate  lines 
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Figure  74  -  Constant  Coordinate  lines 


63 


PHYSICAL  SPACE  GRID 

Ml-0.70,  RE1-200K/FT,  ETR-O.OS,  GRIO:  128  X  32 


Figure  75  -  Constant  Coordinate  lines 
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Figure  76  -  Full  Field  Pressure  Contours 


Figure  77  -  Full  Field  Velocity  Vector  Contours 
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Figure  78  -  Full  Field  Unit  Mass  Flow  Contours 
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Figure  79  -  Full  Field  Mach  Number  Contours 
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Figure  80  -  Pressure  Contours 


Figure  81  -  Velocity  Vector  Contours 
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Figure  82  -  Unit  Mass  Flow  Contours 


Figure  83  -  Mach  Number  Contours 
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Figure  84  -  Converged  Euler  Solution 
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Figure  85  -  Surface  Pressure  Distribution 
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Figure  86  -  Surface  Mach  Number  Distribution 
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SURFACE  OENSITY  DISTRIBUTION 
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Figure  87  -  Surface  Density  Distribution 
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Figure  88  -  Surface  Pressure  Coefficient  Distribution 
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Figure  89  -  Surface  Mach  Number  Distribution 
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Figure  90  -  Surface  Density  Distribution 
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Figure  91  -  Full  Field  Velocity  Profiles 
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Figure  93  -  Full  Field  Unit  Mass  Flow  Profiles 
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Figure  95  -  Full  Field  Pressure  Contours 
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Figure  96  -  Full  Field  Velocity  Contours 
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Figure  97  -  Full  Field  Unit  Mass  Flow  Contours 
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Figure  98  -  Full  Field  Mach  Number  Contours 
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Figure  99  -  Pressure  Contours 


Figure  100  -  Velocity  Vector  Contours 
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Figure  101  -  Unit  Mass  Flow  Contours 
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Figure  102  -  Mach  Number  Contours 
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Figure  103  -  Surface  Pressure  Coefficient  Distribution 
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Figure  104  -  Surface  Mach  Number  Distribution 
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SURFACE  DENSITY  DISTRIBUTION 
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Figure  105  -  Surface  Density  Distribution 
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Figure  106  -  Surface  Pressure  Coefficient  Distribution 
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Figure  107  -  Surface  Mach  Number  Distribution 
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Figure  108  -  Surface  Density  Distribution 
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Figure  109  -  Full  Field  Density  Profiles 
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Figure  111  -  Full  Field  Energy  Profiles 
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Figure  112  -  Full  Field  Unit  Mass  Flow  Profiles 
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Figure  114  -  Full  Field  Pressure  Contours 


Figure  115  -  Full  Field  Velocit  Vector  Contours 
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Figure  116  -  Full  Field  Unit  Mass  Flow  Contours 
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Figure  117  -  Full  Field  Mach  Number  Contours 
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Figure  118  -  Pressure  Contours 


Figure  119  -  Velocity  Vector  Contours 
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Figure  120  -  Unit  Mass  Flow  Contours 


Figure  121  -  Mach  Number  Contours 
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Figure  122  -  Surface  Pressure  Coefficient  Distribution 


SURFBCE  MflCH  NO.  DISTRIBUTION 

MI-0. 81,  REWOOK/fT,  ETfl-0.05,  GRIO:  128  X  32 


oj 


Figure  123  -  Surface  Mach  Number  Distribution 


SURFACE  DENSITY  DISTRIBUTION 

Ml-0.81,  RE1-200K/FT,  ETR-0.05,  GRIO:  128  X  32 


Figure  124  -  Surface  Density  Distribution 


SURFACE  CP  DISTRIBUTION 

Ml -0.81,  REI-200K/FT,  GRIO:  128  X  32,  ITLOCL-3 


Figure  125  -  Surface  Pressure  Coefficient  Distribution 
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Figure  126  -  Surface  Mach  Number  Distribution 


SURFACE  DENSITY  DISTRIBUTION 
MI-0.8'1,  REl-200K/nr,  eTfl-O.OS,  GRID:  128  X 


32 


Figure  127  -  Surface  Density  Distribution 
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Figure  129  -  Full  Field  Velocity  Profiles 
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Figure  131  -  Full  Field  Unit  Mass  Flow  Profiles 
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LOG  10  RHO  RESIOUflL  I TERRT I ON  HISTORY 
Hl-0.84,  RE1-200K/FT,  ETA-O. OS,  GR10:  128  X  32 


Figure  132  -  Residual  Iteration  Time  History 
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SURFACE  Cp  DISTRIBUTION 

ITER=4K,  CN=0.50,  Ml=0.84,  ETA=0.05 
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Figure  133  -  Euler  Solution  Surface  Pressure  Distribution 


PRESSURE  CONTOURS 
ITER=4K,  CN=0.50,  Ml=0.84,  ETA=0.05 


X 

Figure  134  -  Euler  Solution  Pressure  Contours 
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